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Abstract. In this paper we propose a new Bayesian estimation method to solve linear 
inverse problems in signal and image restoration and reconstruction problems which has 
the property to be scale invariant. In general, Bayesian estimators are nonlinear functions 
of the observed data. The only exception is the Gaussian case. When dealing with linear 
inverse problems the linearity is sometimes a too strong property, while scale invariance 
often remains a desirable property. As everybody knows one of the main difficulties 
with using the Bayesian approach in real applications is the assignment of the direct 
(prior) probability laws before applying the Bayes' rule. We discuss here how to choose 
prior laws to obtain scale invariant Bayesian estimators. In this paper we discuss and 
propose a familly of generalized exponential probability distributions functions for the 
direct probabilities (the prior p(x) and the likelihood p(y\x)), for which the posterior 
p(x\y), and, consequently, the main posterior estimators are scale invariant. Among 
many properties, generalized exponential can be considered as the maximum entropy 
probability distributions subject to the knowledge of a finite set of expectation values of 
some knwon functions. 



1. Introduction 

We address a class of linear inverse problems arising in signal and image recon- 
struction and restoration problems which is to solve integral equations of the form: 



9ij 



f /(rOMOdr' + fetf, i,j = l, — ,M, (1) 
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where r' G R 2 , f(r') is the object (image reconstruction problems) or the original 
image (image restoration problems), gij are the measured data (the projections in 
image reconstruction or the degraded image in image restoration problems), bij 
are the measurement noise samples and hij(r') are known functions which depend 
only on the measurement system. To show the generality of this relation, we give 
in the following some applications we are interested in: 
— Image restoration: 

ff i = 1 • • • N 

9(x2,yj) = Jj f(x ,y)h(xi-x ,yj-y)dx dy +b(xi,yj) , j = {...' M > 



where g(xi, yj) are the observed degraded image pixels and h(x, y) is the point 
spread function (PSF) of the measurement system. 



— X-ray computed tomography (CT): 



ff i = 1 • • • N 

g(ri,<f>j)= f(x,y)6(r i -xcos^ i -ysm(l)i)dxdy+b(r i ,(() j ) , - = 1 '...' M 

where g(n,(j>j) are the projections along the axis ri = x cos (f>i — ysin 
having the angle <pj , and which can be considered as the samples of the Radon 
transform (RT) of the object function f(x,y). 

Fourier Synthesis in radio astronomy, in SAR imaging and in diffracted wave 
tomographic imaging systems: 



9(uj,Vj) = JJ f(x,y) exp [-j(ujX + v j y)] dxdy + b(u j ,v j ), j = 1, ■ • -,M, 

where Uj = (uj,Vj) is a radial direction and g(tij,Vj) are the samples of the 
complex valued visibility function of the sky in radio astronomy or the Fourier 
transform of the measured signal in SAR imaging. 

Other examples can be found in ||, [?], [|, || [|. 

In all these applications we have to solve the following ill-posed problem: how 
to estimate the function f{x,y) from some finite set of measured data which may 
also be noisy, because there is no experimental measurement device, even the most 
elaborate, which could be entirely free from uncertainty, the simplest example 
being the finite precision of the measurements. 

The numerical solution of these equations needs a discretization procedure 
which can be done by a quadrature method. The linear system of equations 
resulting from the discretization of an ill-posed problem is, in general, very ill- 
conditioned if not singular. So the problem is to find a unique and stable solution 
for this linear system. The general methods which permit us to find a unique and 
stable solution to an ill-posed problem by introducing an a priori information on 
the solution are called regularization . The a priori information can be either in 
a deterministic form (positivity) or in a stochastic form (some constraints on the 
probability density functions). 

When discretized, these problems can be described by the following: 

"Estimate a vector of the parameters x G R™ (pixel intensities in an im- 
age for example) given a vector of measurements y € R m (representing, 
for example, either a degraded image pixel values in restoration prob- 
lems or the projections values in reconstruction problems) and a linear 
transformation A relating them by: 

y = Ax + b, (2) 

where b represents the discretization errors and the measurement noise 
which is supposed to be zero-mean and additive." 

In this paper we propose to use the Bayesian approach to find a regularized so- 
lution to this problem. Noting that the Bayesian theory only gives us a framework 
for the formulation of the inverse problem, not a solution of it. The main difficulty 



is, in general, before the application of the Bayes' formula, i.e.; how to formulate 
appropriately the problem and how to assign the direct probabilities. Keeping 
this fact in mind, we propose the following organization to this paper: In section 
2. we give a brief description of the Bayesian approach with detail calculations of 
the solution in the special case of Gaussian laws. In section 3. we discuss about 
the scale invariance property and propose a familly of prior probability density 
functions (pdf) which insure this property for the solution. Finally, in section 4., 
we present some special cases and give detailed calculations for the solution. 



2. General Bayesian approach 

A general Bayesian approach involves the following steps: 

— Assign a prior probability law p(x) to the unknown parameter to translate 
our incomplete a priori information (prior beliefs) about these parameters; 

— Assign a direct probability law to the measured data p(y\x) to translate the 
lack of total precision and the inevitable existence of the measurement noise; 

— Use the Bayes' rule to calculate the posterior law p(x\y) of the unknown 
parameters; 

— Define a decision rule to give values x to these parameters. 

To illustrate the whole procedure, let us to consider an example; the Gaussian case. 
If we suppose that what we know about the unknown input x is its mean E {x} = 
Xq and its covariance matrix E{(:r — Xq){x — Xo)*} = R x = cr^P, and what we 
know about the measurement noise b is also its covariance matrix E { bi>* } = Rb = 
a? I, then we can use the maximum entropy principle to assign: 



p(x) oc exp 

and 



■^(y ~ Ax) t R b ~ 1 {y - Ax) 



(3) 
(4) 



p(y\x) oc exp 

Now we can use the Bayes' rule to find: 

p(x\y) ocp(y\x)p(x), (5) 

and use, for example, the maximum a posteriori (MAP) estimation rule to give a 
solution to the problem, i.e.; 

x = argmax{p(£c|y)} , (6) 

Other estimators are possible. In fact, all we want to know is resumed in the 
posterior law. In general, one can construct a bayesian estimator by defining a 
cost (or utility) function C(x, x) and by minimizing its mean value 

x = argrmn {E X |y {C(z, x)}} = argrmn | J C(z, x)p(x\y) dcc| . 

The two classical estimators: 



— Posterior mean (PM): x = Ex\y {x} = J xp(x\y) das, 

is obtained when defining C(x, x) = (x — x) l (x — x), and 

— Maximum a posteriori (MAP): x = argrnax {p(x\y)}, 

is obtained when defining C(x, x) = 1 — S(x — x). 

Now, let us go a little further inside the calculations. Replacing (||), and (|J) 
in (||), we calculate the posterior law: 

with J(x) — (y—Ax) t (y-Ax)+\(x-x ) t P~ 1 (y~xo), 

where A = a 2 /a 2 .. The posterior is then also a Gaussian. We can know use any 
decision rule to obtain a solution. For example the maximum a posteriori (MAP) 
solution is obtained by: 

x = argmax{p(x|y)} = argmin{J(a:)} . (7) 

Note that in this special Gaussian case both estimators, i.e.; the posterior mean 
(PM) and the MAP estimators are the same: 

x = E x \y {x} = argrnax {p(:r|y)} (8) 
and the minimization of the criterion J(x) which can also be written in the form: 
J(x) = \\y - Ax\\ 2 + X\\x - x Q \\ 2 p (9) 

can be considered as a regularization procedure to the inverse problem (||). Indeed, 
the Bayesian approach will give us here a new interpretation of the regularization 
parameter in terms of the signal to noise ratio, i.e.: A = cr^/cr^. 

J(x) is a quadratic function of x. The solution x is then a linear function of the 
data y. This is due to the fact that the problem is linear and all the probability 
laws are Gaussian. Excepted this case, in general, the Bayesian estimators are 
not linear functions of the observations y. However, we may not need that the 
solution be a linear function of the data y, but the scale invariance is the minimum 
property which is often needed. 



p(x\y) oc exp 
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3. Scale invariant Bayesian estimators 

What we are proposing in this paper is to study in what conditions we can obtain 
estimators who are scale invariant. Note that linearity is the combination of 

additivity: < ^ 1 =>■ y 1 + y 2 i— > X\ + X2, 

[ y-2 l— * x 2 

and 

scale invariance: y 1 <— * X\ => Vfc > 0, ky 1 t— > kx\. 



In a linear inverse problem what is often necessary is that the solution be scale 
invariant. As we have seen in the last section when all the probability laws are 
Gaussian then the Bayesian estimators are linear functions of the data, so that 
the methods based on this assumption have not to take care about the scale of the 
measured data. The Gaussian assumption is very restrictive. On the other hand, 
more general priors yield the Bayesian estimators which are nonlinear functions of 
data, so the result of the inversion method depend on the absolute values of the 
measured data. In other words, two users of the method using two different scale 
factors would not get the same results, even rescaled: 
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A general nonlinear (scale variant) estimation method 



What we want to specify in this paper is a family of probability laws for which these 
estimators are scale invariant. So the user of the inversion method can process the 
data without worrying about rescaling them to an arbitrary level and two users of 
the method at two different scales will obtain the proportional results: 
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A scale invariant estimation method 



To do this let us note 

— 6 all the unknown parameters defining our measuring system (noise variance 
a 2 and the prior law parameters for example), 

— p 1 (xi\y 1 ; 61) and p k {x k \y k ;9k) the two expressions of the posterior law for 
scale 1 and for scale k with 

x k = feaji, y k = ky 1 . 

Then, what we need is the following: 



39 k = f{6i,k)\Vk>0,Vx 1 ,y 1 , p k {x k \y k ; B k ) = — pi{x!\ Vl ; 0i), (10) 



which means that the functional form of the posterior law remains unchanged 
when the measurement's scale is changed. Only we have to modify the parameters 
k = f(0\,k) which is only a function of B\ and the scale factor k. 



However, not all estimators based on this posterior will be scale invariant. The 
cost function must also have some property to obtain a scale invariant estimator. 
So, the main result of this paper can be resumed in the following theorem: 

Theorem: If 30 k = f(6 1 ,k)\Vk > 0,Vx 1 ,y 1 , 

Pk (x k \y k ;8 k ) = ^p 1 {x 1 \y l ;6 1 ), 

then any bayesian estimator with a cost function C(x, x) satisfying: 

C(x k ,x k ) = a k + b k C(x,x), 
is a scale invariant estimator, i.e.; 

x k (y k ;0 k ) = kxiiy^O^. 

Proof: In fact, it is easy to see the following: 

x k (y k ;O k ) = argrnm j J C(z k ,x k )p k (x k \y k ;O k ) dx fe | 

= k argmin j f[b k C(z 1 ,x 1 ) + a k }^- p 1 (x 1 \y 1 ; 0i)fc" dxj 
z-t \J k n 

= k argmin \b k J C(z 1 ,x 1 )p 1 (x 1 \y 1 ;0 1 )dx 1 +a fc | 

= k argmin <^J C(z 1 ,x 1 )p 1 (x 1 \y 1 ;0 1 ) dx^j 
= kx^y^Oi) 

Note the great significance of this result, even if the estimateur x(y;0) is a 
nonlinear function of the observations y it stays scale invariant. 

Now, the task is to search for a large familly of probability laws p(x) and p(y\x) 
in a manner that the posterior law p{x\y) remains scale invariant. We propose to 
do this search in the generalized exponential familly for two reasons: 

— First the generalized exponential probability density functions form a very 
rich one, and 

— Second, they can be considered as the maximum entropy prior laws subject 
to a finite number of constraints (linear or nonlinear). 

Noting also that if p(x) and p(y\x) are scale invariant then the posterior p(x\y) 
is also scale invariant and that there is a symmetry for p(x) and p(y\x), so that 
it is only necessary to find the scale invariance conditions for one of them. In the 
following, without loss of generality, we consider the case wherep(y|a;) is Gaussian: 

p(y\x;a 2 ) cx cxp [-x 2 {x, y: a 2 )] , with x 2 (x,y;<J 2 ) = iv ~ Hx fly - Hx }> 

(11) 



and find the conditions for p(x) to be scale invariant. We choose the generalized 
exponential pdf's for p(x), i.e.; 



p(x; A) cx exp 



\j(j)i{x) 



(12) 



and find the conditions on the functions <fii(x) for which p(x) is scale invariant. 

Note that these laws can be considered as the maximum entropy prior laws if 
our prior knowledge is: 

— What we know about x is: 



E{<f>i(x)} = d h 
and what we know about the noise b is: 



1 



5 ) ' J 



E{6} = 0, 

E{6b'} =R b = a 2 I, 

where R b is the covariance matrix of b. 

Now, using the equations ( pi] ) and ( |l2| ) and noting by 6 = (er 2 , Ai, 
A = (Ai, • • • , A r ), and by (f>{x) — (cj) 1 (x), • ■ • , <f) r (x)), we have 

p(x\y; 0) cx exp [~x 2 (a;, y; a 2 ) - A*0(«)] > 

and the scale invariance condition becomes: 

Vfc > O.VaJi.tfx, Xfc(a5fc.y*;o|)+^(xfc) =x 2 (xi,y 1 ;(rf)+A t 1 0(xi) + cte. 

But with the Gaussian choice for the noise pdf we have 



A r ), by 
(13) 



Vfe > o.Vxi.wx, xt(x k ,y k ;ol) = 7^2 \\yk-Hx k \\ 2 = 

and so the condition becomes 



fc 2 \\ yi -Hxi\\ 2 = xl{xi,y 1 \a 2 1 ), 



Vfc > 0,Vx, \\<p(xk) = X\cj)(xi) + cie, 



(14) 



or equivalently, 



Pk(xk; Afe) = — pi(xi; Ai) with \ k = f(\i,k). 

Thus, in the case of centered Gaussian pdf for the noise, to have a scale invariant 
posterior law it is sufficient to have a scale invariant prior law. 
Now, assuming interchangeable (independent) pixels, i.e.; 



p(x; A) = exp 



A + ^2 ^i<t>i( x ) 



N 



(15) 



or equivalently, 

N 

= J^&fo) (16) 

3=1 

we have to find the conditions on the scalar functions fa (x) of scalar variables x 
who satisfy the equation (|lj) or equivalently 

r r 

Vfc>0,Va;, ^2\i(k)fa(kx)=^2Xi(l)fa(x)+cte (17) 

i=l i=l 

We have shown (see appendix) that, the functions fa(x) which satisfy these con- 
ditions are all either the powers of x or the powers of In x or a multiplication of 
them. The general expressions for these functions are: 

M /N m -1 \ N M 

fax) = I £ c mn (lnx) n j a^+J c „(lnx)", with M<r and ^ 7V m = 

rn — 1 \ n — / n— m— 

(18) 

where M and iV m are integer numbers, and c m „, co n and a m are real numbers. 
For a geometrical interpretation and more details see appendix. The following 
examples show some special and interesting cases. 

One parameter laws: Consider the case of r = 1. In this case we have 

p(x; A) oc exp [— Xfax)] . (19) 

Applying the general rule with 

(M = 0,N = 1, — >c 00 + c ilnx 
r \ M = 1,N = 0, N 1 = l, — » coo + ci a; Ql 

we find that the only functions who satisfy these conditions are: 

fa x )\ = |a; Q ,lna;| (20) 

where a is a real number. There isc two interesting special cases: 

— fax) = x a , resulting to: p{x) oc exp [— Ax"] , a > 0, A > 0, which is a gener- 
alized Gaussian pdf, and 

— fax) — In a;, resulting to: p(x) oc exp [— A In a;] , which is a special case of the 
Beta pdf. 

Note that the famous entropic prior law: p{x) oc exp [—Ace In x] of Gull and 
Skilling jll], U docs not verify the scale invariance property. But, if we add one 
more parameter 

p(x) oc exp [—Ax In a; + fix] , 
then, it will satisfy this condition as we can see in the next section. 



Two parameters laws: This is the case where r = 2 and we have: 

p(x; A) oc cxp [— A0i(x) - /U0 2 (x)] , (21) 
and applying the general rule: 



M = 2, N 


= 0,N, 


= 1,JV 2 = 1, - 


— » c 00 - 


f- CioX" 1 - 


f- c 20 x« 2 


M = 1,N 


= 0,JVi 


= 2, 


— > c o - 


I" Ciox"! - 


- c\\x ai In a; 


M = 1,N Q 


- l.JVi 


= 1, 


— » coo - 


h CloX" 1 - 


1- coi In a; 


M = 0, N 


-2, 




— » coo - 


1- coi In a; 


+ c 2 In 2 x 



we see that in this case the only functions (0i, 02 ) which satisfy these conditions 
are: 



| (^1 (z) , 02 (ar)) | = | (x ai , a;" 2 ) , (x" 1 , x Ql In x) , (x" 1 , In x), (In x, In 2 x) j 

(22) 

where ai and a 2 are two real numbers. Special cases are obtained when we choose 
02 (x) = x, the only possible functions for 0i(x) are then: 

{x a ,mx,xlnx}. (23) 

and we have the following interesting cases: 

— 01 (x) = x 2 , resulting to: p(x) oc exp [—Ax 2 — fix] oc exp —A (x + ^j) 2 , 

which is a Gaussian pdf W (m = cr 2 = ^) . 

— 0i (x) = In x, resulting to: p(x) oc exp [— A lnx — fix] = x~ A cxp [— fix] , which 
is the Gamma pdf, and finally, 

— 0i (x) = xlnx, resulting to: p(x) oc cxp [—Ax lnx — fix] . which is known as 
the entropic pdf. 

Three parameters laws: This is the case where r = 3. Once more applying the 



r = 3 



which means: 

(0i(x),0 2 (x),0 3 (x)) 



we 


find: 












M 


= 3,7V 


= O.JVi 


= 1,^2 


= 1,^3 = 1, 


-> c 00 - 


h Ci X Ql + c 20 x" 2 + c 30 x a3 


M 


= 2,7V 


= 0,iVi 


= 1,JV 2 


-2, 


—> coo - 


h ci x ai + c 20 x a2 + c 2 ix" 2 lnx 


M 


= 2,7V 


- l.JVi 


- 1,JV 2 


= 1, 


-» c 00 - 


- coi lnx + ci x Ql + c 2 ox Q2 


M 


= l,JVo 


= 0,JVi 


-3, 




-» c 00 - 


- ciqx" 1 + cnx Ql lnx + c i2 x ai In 2 x 


M 


= 1,JV 


= 1, JVi 


= 2, 




— > cqo + c oi lnx + ciqx" 1 + cux ai lnx 


M 


= 1,JV 


= 2,JVi 


= 1, 




-> c 00 - 


- coi In x + c 02 In 2 x + ciox" 1 


M 


= 0,7V 


= 3, 






-> c 00 - 


f- c i In x + c 02 In 2 x + c 3 In 3 x 



(x Ql ,x Q2 ,x" 3 ), (x ai ,x Q2 ,lnx), (x ai ,x Ql lnx,x Ql ln 2 x), 
(x Ql , x" 1 lnx, lnx), (x" 1 , x" 2 , x" 2 lnx), (x Ql , lnx, In 2 x), 

(lnx, In 2 x, In 3 x) 
(24) 



where ai, a 2 and 0:3 are three real numbers. 



4. Proposed method 

The general procedure of the inversion method we propose can be resumed as 
follows: 

— Choose a set of functions fa (x) between the possibles ones described in the last 
section and assign the prior p(x). In many imaging applications we proposed 
and used successfully the following two parameters one: 

N 

p(x;X) oc exp [— XiH(x) — \2S(x)] , with H(x) = \^(f>i(xj), and S(x) - 

3=1 

where <j>\ (x) and (f>2 (x) choosed between the possible ones in (^2|) or (^) . 

— When what we know about the noise b is only its covariance matrix E {fob*} 
Rb = <7? I, then using the maximum entropy principle we have: 

with Q(x) = (y - Ax) t R b ~ 1 {y - Ax). 

We may note that p(y\x) is also a scale invariant probability law. 

— Using the Bayes' rule and MAP estimator the solution is determined by 

x = argm&x. {p(x\y)} = argmin {J(x)} , with J(x) = Q(x)+\iH(x)+\2S(x). 

Note here also that, for the cases where one of the functions <j>i(x) or <j>2 (x) 
is a logarithmic function of x, we have to constraint its range to the positive 
real axis, and we have to solve the following optimization problem 

x = argmax{p(a;|y)} = argmin {J(x)} . 

This optimization is achieved by a modified conjugate gradients method. 

— The choice of the functions <fti{x) and the determination of the parameters 
(Ai, A2) in the first step is still an open problem. 

In imaging applications we propose to do this choice from our prior knowledge 
on the nature of interested quantity (physics of the application). For example, 
if the object a; is a real quantity equally distributed on the positive and the 
negative reals then a Gaussian prior, i.e.; (<fii(x) = x,4>2{x) — x 2 ) is conve- 
nient. But, if the object a; is a positive quantity or if we know that it represents 
small extent, bright and sharp objects on a nearly black background (images in 
radio astronomy, for example), then we may choose ((f>i(x) = x, 4>2(x) = lnx), 
or (0i (x) = x, ^{x) = xlnx) which are the priors with longer tails than the 
Gaussian or truncated Gaussian one. 

When the choice of the functions (0i(x), <p2 (x)) is done, we still have to deter- 
mine the hyperparameters (Ai, A2). For this two main approaches have been 
proposed. The first is based on the generalized maximum likelihood (GML) 
which tries to estimate simultaneously the parameters x and the hyperparam- 
eters 9 — (Ai, A2) by 

(x, 9) = arg max {p(x, y; 9)} = arg max {p(y\x) p(x; 9)} , (25) 
(x,9) (x,0) 



N 



p(y\x) oc exp 



-M x ) 



and the second is based on the marginalization (MML), in which the hyper- 
parameters 9 are estimated first by 

6 = argmax 0) — j p(x, y; 6) dx X = argmax^ J p(y\x)p(x; 0) Ax X , 

° ' (26) 

and then used for the estimation of x: 

x = argmax |p(a;|y; <?) j = argmax |p(y|a;)p(a;|0) j . (27) 

What is important here is that both methods preserve the scale invariant 
property. For practical applications we have recently proposed and used a 
method based on the generalized maximum likelihood J|, [| which has been 
successfully used in many signal and image reconstruction and restoration 
problems as we mentionned in the introduction [1(J . 

5. Conclusions 

Excepted the Gaussian case where all the Bayesian estimators are linear functions 
of the observed data, in general, the Bayesian estimators are nonlinear functions 
of the data. When dealing with linear inverse problems linearity is sometimes a 
too strong property, while scale invariance often remains a desirable property. In 
this paper we discussed and proposed a familly of generalized exponential proba- 
bility distributions for the direct probabilities (the prior p(x) and the likelihood 
p(y\x)), for which the posterior p(x\y), and, consequently, the main posterior es- 
timators are scale invariant. Among many properties, generalized exponential can 
be considered as the maximum entropy probability distributions subject to the 
knowledge of a finite set of expectation values of some knwon functions. 

1. Appendix: General case 

We want to find the solutions of the following equation: 

r r 

Vfc>0,Vx, Y,\ i (k)<p i {kx)=J2\ i (l)(P i (x)+P(k) (A.l) 

i=l i=l 

Making the following changes of variables and notations 

1/k = k, kx = x, Xi(k) = A^(fc), and Pi{k) — f3i(k), 



equation (A.l) becomes 

r r 

A, (fc)& (x) = J2 Ml)* (kx) + P(k) 

For convenience sake, we will drop out the tilde ~, and note A^(l) = A 2 ;, so that we 
can write 

r r 

5^Ai(k)&(aO=5^Ai(fe(fcaO + /9(k) 



Noting 



we have 



S(x) — ^""^ Xj<f>j(x), and so S(kx) = Xj(j>i(kx) 



5^Ai(fc)^(i) = S(fcr) + 



Deriving r — 1 times this equation with respect to k we obtain 



i=l 
r 

x;Aj'(fc)^(x) 



a;S"(fca;) +/3'(fc) 
a; 2 5"(fcx)+/3"(A;) 



(A.2) 



(A.3) 



- r - 1 5( r - 1 )(fet)+^ f - 1 )(fc) 



Combinig equations (|A.2|) and (A.3) in matrix form we have 



/ A x (fc) 
A'/(fe) 



A r (fe) \ 
A^(fc) 

A;'(fc) 



02 (a:) 

03 (a;) 



/ S(fcc) + /3(A) \ 

xS"(fex) + /3'(/c) 
x 2 ^"(fca;)+/3"(A:) 



4 r-1) (*)/ \^.(as)/ \a; r - 1 1 5('"- 1 )(A;ar)+ ( 9('- 1 )(A;)/ 

(A.4) 

If this matrix equation can be inverted, this means that any function 0j(x) is 
a linear combination of S(kx) + f3(k) and its (r — 1) derivatives with respect to k: 



fcix) =J2m(k) Ix^S^ikx) + fli-VQe) 



i=0 



(A.5) 



and if this is not the case, this means that there exists an interval for k, for which 
some of the functions Xi(k) are linear combinations of the others [0. In this case 
let us show that we will go back to the situation of the problem of lower order r. 
Let us to assume that the last column of the matrix is a linear combination of the 
others, i.e.; 



X r (k) = J2 7i Xi(k). 



Putting this in the equation (A.l) will give 



r-l 



Xi(k)4>i(kx)- 



y^7iAi(fc) 



r-l 



.{kx)= Y"A i (l)0 i (a ; )+^(fc)- 



4> r (x) 



and noting ipi(x) — 4>i(x) + ji4> r (x) and ipi(kx) — 4n{kx) + r )i<j) r {kx) we obtain 



{k)1>i (kx) =j2x i (l)1> i (x)+ p(k) 



i=l 



which is an equation in the same form of (A..1), but of lower order. 

Deriving now both parts of the equation ( |A.5| ) with respect to k and noting 
kx — u we obtain 



y j at u l S' l (u) = a 



(A.6) 



i=0 



This is the general expression of a rth order Euler-Cauchy differential equation 
H, H which is classically solved through the change of variable u — e x , and one 
can find the general expression of its solution in the following form: 



M /JVm-l 



No 



S(x) = [ cmn(lna;) n c „(lnx)" with M = 0, • • -r, and N " 



m— 1 \ n— 



71=0 



m=0 



(A.7) 

where M and N m are integer numbers, and c mn , co n and a m are real numbers. In 
fact the most general solution also incorporate terms of the form 



yj(lnx)™ (a n cos(lna;) + (3 n sin(ln.T)) 



derived from complex a m and c mn . But we will not consider these terms because 
the resulting pdf's have oscillatory behavior around zero. 



One can give a geometric interpretation of the solutions given in (A.7). For 
any given order r make a (r + 1) x (r + 1) table in the form 



In' x 


















ln z x 










In a; 












1 


X 












1 


x ai 









and let r mass points fall down into the columns. To each filled box is assigned a 
function </>j (x) by multiplying the corresponding powers of x and In x on the same 
line and the same column. To illustrate this, we give in the following the three 



first cases: 



Case r = 1: 



In a; 


b 




1 


X 


a 




1 





Case r = 2: 



h\ z x 


d 






In a; 


bd 


c 




1 


X 


abc 


a 




1 


x ai 


x a 2 



Case r = 3: 



\n 6 x 


9 








In'x 


fg 


c 




In a; 


bdfg 


dc 


e 




1 


X 


abcdef 


abe 


a 




1 


x ai 


x a2 


a;" 3 




c 


>i(z) 


<h(%) 


<h( 


r) 


a 




x a i 


x a2 


/ 3 


b 




x ai 


x a2 


In a; 


c 




x ai 


x ai In x 


a;" 1 In 2 a; 


d 




x ai 


x ai In a; 


In a; 


e 




x ai 


a;" 2 


x a2 In a; 


f 




x ai 


In a; 


In 2 


X 


9 




In a; 


In 2 a; 


In 3 







cj>{x) 


a 


x a ^ 


b 


In a; 



01 (x) <j) 2 (x) 



In a; 



In a; 
a;" 1 In a; 
In 2 a; 
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